#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <vector>
#include "es-bvs-bl-common.h"

#define DEBUG

//
//  benchmarking program
//
int main( int argc, char **argv )
{    
  if( find_option( argc, argv, "-h" ) >= 0 )
  {
    printf( "Options:\n" );
    printf( "-h to see this help\n" );
    printf( "-n <int> to set the number of particles\n" );
    printf( "-o <filename> to specify the output file name\n" );
    printf( "-px <int> -py <int> to set the number of pthreads\n"); 
    return 0;
  }
  
  int n = read_int( argc, argv, "-n", 1000 );
  int NUM_THREADS_X = read_int( argc, argv, "-px", 2 );
  int NUM_THREADS_Y = read_int( argc, argv, "-py", 2 );
  char *savename = read_string( argc, argv, "-o", NULL );
  
  FILE *fsave = savename ? fopen( savename, "w" ) : NULL;
  particle_t *particles = (particle_t*) malloc( n * sizeof(particle_t) );
  double grid_size;
  set_size( n, grid_size );
  
  // This should probably be parallelized
  init_particles( n, grid_size, particles );

  const double cell_size = MIN_CELL_SIZE;
  const int bins_dim = (int)((grid_size+cell_size)/cell_size);
  const int bins_thread_size_x = bins_dim / NUM_THREADS_X 
				  + (bins_dim % NUM_THREADS_X > 0 ? 1 : 0);
  const int bins_thread_size_y = bins_dim / NUM_THREADS_Y 
				  + (bins_dim % NUM_THREADS_Y > 0 ? 1 : 0);

#ifdef DEBUG
  printf("grid_size = %lf\n",grid_size);
  printf("threadx = %d thready = %d\n",NUM_THREADS_X,NUM_THREADS_Y);
  printf("array of bins is = %d by %d\n",bins_dim, bins_dim);
  printf("each thread manages bins of size %d by %d\n", bins_thread_size_x, bins_thread_size_y);
#endif
  std::vector< particle_t* >* bins[bins_dim*bins_dim];
  for ( int i = 0; i < bins_dim*bins_dim; i++ ){
    bins[i] = new std::vector< particle_t* >();
  }
  
  int interact_neighbors[9] = {(-bins_dim-1), (-bins_dim), (-bins_dim+1),
			       -1, 0, 1,
			       (bins_dim-1), (bins_dim), (bins_dim+1)};
  
  double simulation_time = read_timer( );
  std::vector< particle_t* > *inter_bin;
  std::vector< particle_t* >::iterator main_part, inter_part;

  //
  //  simulate a number of time steps
  //
  #pragma omp parallel num_threads(NUM_THREADS_X*NUM_THREADS_Y)
  for( int step = 0; step < NSTEPS; step++ ) {
    //
    //  bin particles
    //
    int bin_x, bin_y, bin_num, thread_x, thread_y, i, proc;
    #pragma omp for private(bin_x,bin_y,bin_num,thread_x,thread_y,i,proc)
    for( proc = 0; proc < NUM_THREADS_X*NUM_THREADS_Y; proc++ ) {
      thread_x = proc % NUM_THREADS_X;
      thread_y = proc / NUM_THREADS_X;
  
      /* Browse entire array with each thread but with variable offsets */    
      for( i = n/(NUM_THREADS_X*NUM_THREADS_Y)*proc; i < n; i++ ){
	bin_x = (int)(particles[i].x/cell_size);
	bin_y = (int)(particles[i].y/cell_size);
	if( (bin_x >= thread_x*bins_thread_size_x) && 
	    (bin_x < (thread_x+1)*bins_thread_size_x) && 
	    (bin_y >= thread_y*bins_thread_size_y) && 
	    (bin_y < (thread_y+1)*bins_thread_size_y) ) {
	  bin_num = bin_y*bins_dim+bin_x;
	  bins[bin_num]->push_back(&particles[i]);
	}
      }
      for( i = 0; i < n/(NUM_THREADS_X*NUM_THREADS_Y)*proc; i++ ) {
	bin_x = (int)(particles[i].x/cell_size);
	bin_y = (int)(particles[i].y/cell_size);
	if( (bin_x >= thread_x*bins_thread_size_x) && 
	    (bin_x < (thread_x+1)*bins_thread_size_x) && 
	    (bin_y >= thread_y*bins_thread_size_y) && 
	    (bin_y < (thread_y+1)*bins_thread_size_y) ) {
	  bin_num = bin_y*bins_dim+bin_x;
	  bins[bin_num]->push_back(&particles[i]);
	}
      }

    }
      
    //
    //  compute forces
    //
    int x,y,bin_to;
    #pragma omp for private(thread_x,thread_y,x,y,i,main_part,inter_bin,bin_to,inter_part)
    for( int proc = 0; proc < NUM_THREADS_X*NUM_THREADS_Y; proc++ ) {
      thread_x = proc % NUM_THREADS_X;
      thread_y = proc / NUM_THREADS_X;
      for( x = thread_x*bins_thread_size_x; 
	   x < min((thread_x+1)*bins_thread_size_x,bins_dim); 
	   x++ ){
	for( y = thread_y*bins_thread_size_y; 
	     y < min((thread_y+1)*bins_thread_size_y,bins_dim); 
	     y++ ){ 
	  i = y*bins_dim+x;
	  for ( main_part=bins[i]->begin(); 
		main_part<bins[i]->end(); 
		main_part++ ) {
	    (*main_part)->ax = 0;
	    (*main_part)->ay = 0;
	    for ( int bin_to = 0; bin_to < 9; bin_to++ )
	    {
	      // Edgar - if bin is on edge this will attempt redundant bin interactions
	      if (i + interact_neighbors[bin_to] < 0 ||
		  i + interact_neighbors[bin_to] >= bins_dim*bins_dim )
		continue;
	      inter_bin = bins[i + interact_neighbors[bin_to]];
	      for ( inter_part=inter_bin->begin(); 
		    inter_part<inter_bin->end(); 
		    inter_part++ ){
		apply_force( **main_part, **inter_part );
	      }
	    }
	  }
	}
      }
    }

    //
    //  move particles and clean out the bins
    //
    #pragma omp for 
    for( int i = 0; i < n; i++ ) {
      move( particles[i], grid_size );
    }
    
    #pragma omp for 
    for( int i = 0; i < bins_dim*bins_dim; i++ ){
      bins[i]->clear();
    }
    //
    //  save if necessary
    //
    #pragma omp master
    if( fsave && (step%SAVEFREQ) == 0 )
      save( fsave, n, particles, grid_size );
    
  }
  simulation_time = read_timer( ) - simulation_time;
  
  printf( "n = %d, simulation time = %g seconds\n", n, simulation_time );
  
  free( particles );
  if( fsave )
    fclose( fsave );
  
  return 0;
}
